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1.  Introduction.  This  survey  of  selected  computational  aspects  of  linear 
algebra  is  addressed  to  members  of  SIAM  who  are  not  specialists  in  numerical 
analysis.  The  reader  is  assumed  to  have  a  general  familiarity  with  the  algebra 
and  analysis  of  finite  vectors  and  matrices,  including  norms,  and  to  know  the 
Gaussian  elimination  process.  A  completely  adequate  background  is  given  in  the 
first  72  pages  of  Faddeeva  [9].  A  much  more  complete  background  for  practical 
matrix  work  is  found  in  Bellman  [3],  Marcus  and  Mine  [38],  and  Wilkinson  [6l]. 

Far  more  extensive  expositions  of  the  computational  methods  of  linear 
algebra  are  to  be  found  in  Fox  [l4],  Noble  [42],  Householder  [28],  and  Wilkinson  [6l] . 

The  author  gratefully  acknowledges  conversations  with  Gene  H.  Golub, 

Richard  Hamming,  and  William  Kahan,  and  especially  the  opportunity  to  see  a 

draft  of  Kahan  [32].  He  also  acknowledges  substantial  debts  to  Cleve  Moler  for  s 

the  use  of  material  from  Forsythe  and  Moler  [l2j. 

2.  Computational  problems  of  linear  algebra.  The  ordinary  computational 
problems  of  linear  algebra  are  concerned  with  matrices  of  real  numbers.  4 

a.  Let  A  be  an  n- rowed,  n-columned  matrix  of  real  numbers.  Let  b  be 
an  n-rowed  column  vector  of  real  numbers.  The  traditional  linear-equations 
problem  is  to  find  an  n-rowed  column  vector  x  such  that 

A  x  -  b  . 

It  is  normally  assumed  that  A  is  a  nonsingular  matrix,  since  then  and  only  then 
does  a  unique  solution  exist  for  all  b. 

b.  With  the  same  A  as  in  part  a,  another  traditional  problem  is  to  find 
the  inverse  matrix  A-^  . 

c.  Let  A  be  an  n-rowed,  n-columned  matrix  of  real  numbers  which  is 
symmetric.  The  third  traditional  problem  is  to  find  some  or  all  of  the  (necessarily 
real)  eigenvalues  of  A.  Recall  that  an  eigenvalue  of  A  is  a  number  k  for  which 
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there  exists  a  column  vector  u  such  that 


A  u  =  X  y  . 

Such  a  vector  u  is  called  a  (column)  eigenvector  of  A  belonging  to  X, 
and  often  the  computational  problem  includes  finding  a  u  belonging  to  each 
eigenvalue  computed..  There  exist  n  orthonormal  eigenvectors  of  A,  one 
belonging  to  each  eigenvalue  of  A. 

d.  Let  A  be  an  unsymmetric  n-rowed,  n-columned  matrix  of  real  numbers. 
Another  traditional  problem  of  linear  algebra  is  to  find  some  or  all  of  its 
eigenvalues,  and  sometimes  also  its  corresponding  column  eigenvectors  and  row 
eigenvectors.  Recall  that  a  row  eigenvector  belonging  to  X  is  an  n-columned 
row  vector  v  such  that 


v  A  =  X  v  . 

When  A  is  not  symmetric,  the  problem  is  complicated  in  many  ways:  First, 

some  of  the  eigenvalues  X  are  ordinarily  complex  numbers.  Second,  there  may 

not  exist  n  linearly  independent  column  eigenvectors,  and  those  which,  exist  are 

not  usually  orthogonal.,/  Indeed,  they  are  likely  to  be  nearly  linearly,  dependent  and 

the  same  holds  for  the  row  eigenvectors.'  Third,,  if  an  eigenvalue  X  is  a  root  of 

multiplicity  k  >  1  of  the  characteristic  eauation  det(A  -  X  i)  =  0,  then 

there  may  exist  anywhere  from  1  to  k  linearly  independent  column  eigenvectors 

belonging  to  X.  (if  A  were  symmetric,  there  would  always  be  k.)  If  the 

number  is  less  than  k,  it  corresponds  to  one  or  more  nondiagonal  blocks  in 

the  Jordan  canonical  form  of  A,  or  equivalently  to  so-called  nonlinear  elementary 

dj visors  of  A.  Fourth,  multiple  or  nearly  multiple  eigenvalues  of  A  are 

likely  to  be  very  rapidly  changing  functions  of  the  elements  a.  .  of  A,  so 

J 

that  computations  are  at  best  tricky. 
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e.  For  any  column  vector  y,  define  the  p-th  power  norm  of  y  to  be 

1 

(1)  ML  =  (  £  |yilp)p  • 

P  i=l 

Here  p  is  a  real  number  with  1  <  p  <  oo  ,  and  y  ,  . ..,  yn  are  the  components 
of  y  in  a  given  coordinate  system.  We  define  the  maximum  norm  as  the  limiting 
case  p  ->  oo  of  (l) ; 

(2)  IIyIIqq  =  max  |y.|  . 

1<  i<  n 


The  norms  most  used  in  numerical  analysis  are  p  =  1,  2,  oo ,  but  statisticians 
are  now  giving  attention  to  values  of  p  between  1  and  2. 

Let  A  be  an  n-rowed,  k-columned  matrix  of  real  numbers,  and  let  b  be 
an  n-rowed  column  vector.  Given  some  p,  a  more  recent  computational  problem 
is  to  find  a  k-rowed  column  vector  x  such  that 


||  A  x  -  b||p  is  minimized  . 

When  p  =  2,  the  usual  case,  this  is  the  linear  least- squares  problem.  For 
p  =  2  the  unit  sphere  in  the  norm  is  very  smooth,  and  methods  of  analysis  work 
well.  However,  for  p  =  1  or  oo  the  unit  sphere  has  many  corners,  and  methods 
of  minimizing  ||Ax  -  b||^  become  combinatorial  or  discrete. 

f .  For  two  n-rowed  column  vectors  x  and  y,  we  define  x  >  y  to  mean 
that  x.^  >  y^  for  all  components  of  x  and  y. 

Let  A  and  b  be  as  in  part  e  above.  Then  an  important  computational  problem 
is  to  describe  the  set  S  of  k-rowed  column  vectors  x  such  that 


I 


A  x  >  b 


S 


I 
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Sometimes,  as  in  linear  programming  problems,  one  looks  for  vectors  x  in  S 
T 

such  that  c  x  is  a  minimum,  where  c  is  a  given  k-rowed  column  vector. 

So  far'  we  have  spoken  only  of  matrices  of  real  numbers.  Similar 
problems  are  posed  occasionally  for  matrices  of  complex  numbers.  Many  of  the 
problems  can  also  be  phrased  for  matrices  whose  elements  are  expressions  in 
indeterminate s  or  letters.  As  methods  of  symbol  manipulation  on  digital  computers 
become  more  accessible  to  computer  users,  problems  of  linear  algebra  with 
matrices  of  letters  will  be  studied  more.  Practical  symbol  manipulation  will 
probably  do  more  to  interest  mathematicians  in  computing  than  anything  that 
has  happened  in  the  computer  era  to  date. 

The  present  discussion  is  limited  to  matrices  of  numbers,  and  moreover 
to  problems  a,  b,  c,  d.  For  discussions  of  problem  e  with  p  =  2,  the 
reader  Is  referred  to  Golub  and  Kahan  [l8l.  For  problem  f  see  presentations 
on  linear  programming  like  Dantzig  (5l. 

Why  do  the  linear  problems  a,  b,  c,  and  d  arise  so  often?  Why  are 
they  important?  The  answer  is  that  linear  operators  are  the  simplest  ones  in 
mathematics,  and  the  only  operators  that  axe  fully  understood  in  principle, 
rience  they  are  a  natural  model  for  an  applied  mathematician  to  use  in  attacking 
■j,  problem.  Even  though  linear  operators  in  infinite-dimensional  spaces  will 
occur  in  analysis  of  differential  equations  (for  example),  the  realities  of 
computing  mean  that  only  finite-dimensional  spaces  can  be  handled  with  digital 
computers. 

More  realistic  models  of  applied  mathematics  are  usually  nonlinear.  But, 
whenever  nonlinear  operators  are  used,  the  actual  solution  of  functional 
equations  almost  always  involves  the  approximation  of  nonlinear  operators  by 
Unear  ones.  A  typical  example  of  this  is  the  use  of  Newton's  method  for  solving 
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a  system  of  nonlinear  equations,  in  which  at  every  step  a  locally  best-fitting 
linear  equation  system  must  be  solved.  Nonlinear  problems  usually  are  very  hard 
In  attacking  them  by  linear  methods,  it  is  essential  th~t  our  linear  tools  be 
very  sharp,  so  th  t  they  can  be  relied  upon  to  work  without  failure.  Only  in 
this  way  can  the  analyst  concentrate  on  the  real  difficulties  of  the  nonlinear 
world.  This  point  of  view  not  only  emphasizes  the  importance  of  being  able  to 
solve  linear  problems,  but  also  the  necessity  of  solving  linear  systems  with 
extremely  reliable  methods. 

Linear  equation  systems  a  arise  directly  mainly  from  two  sources.  One 
is  from  an  approximation  to  linear  functional  equations,  usually  ordinary  or 
partial  differential  equations.  The  other  source  is  a  problem  of  data  fitting, 
interpolation,  or  approximation  by  linear  families  of  functions. 

Eigenvalue  problems  usually  arise  from  studies  of  vibration  or  stability 
or  resonance  of  linear  physical  systems  (e.g.,  .flutter  of  aircraft  and  criti¬ 
cality  of  reactors),  or  from  factor  analysis  problems. 

An  excellent  textbook  by  Noble  [42]  gives  a  number  of  physical  examples 
of  computational  matrix  problems. 


3.  A  closer  look  at  the  problems.  Since  actual  computers  have  finite 

storage  capacity  and  a  finite  precision,  we  need  to  have  a  closer  look  at  the 

nature  of  the  matrices  A  and  the  computational  problems. 

Is  the  matrix  A  dense  (most  elements  a^  ^  ^  0),  or  is  it  sparse 

(most  elements  a.  .  =  0)?  If  A  is  sparse,  do  the  nonzero  elements  form  a 

i>  J 

significant  pattern?  For  example,  is  A  triangular  (a..  ^  =  0  for  i  >  j  or 

for  i  <  j)?  Is  it  of  Hessenberg  form  (ai  ^  =  0  for  i  >  j  +  1  or  for 

j  >  i  +  l)?  Is  it  a  band  matrix  (a^  j  --  0  for  |i  -  j  |  >  m,  where  m  «  n? 

Is  it  a  tridiagonal  matrix  (i.e.,  a  band  matrix  with  m  --  l)?  All  these  special 
forms  occur  frequently,  and  can  be  given  special  consideration. 

5 


Is  the  matrix  A  symmetric?  Positive  tie  f  ini  to?  If  it  is  sparse,  is  ti.e 

pattern  associated  with  the  adjacency  matrix  of  sane  graph?  Frequently  matricc 
associated  with  structures  or  with  partial  difference  equations  are  best  under¬ 
stood  in  terms  of  the  associated  graph. 

Are  the  elements  aJ  stored  in  the  computer  memory,  to  be  retrieved 

i>  d 

when  needed,  or  are  they  regenerated  from  some  algorithm,  as  needed?  One  might 
define  the  informational  content  of  a  matrix  as  the  number  of  cells  needed 
(on  a  certain  computer)  to  store  the  data  and  program  to  obtain  all  the  a.  . . 
The  author  knows  of  no  work  on  this  concept,  which  is  clearly  relevant  to 
matrix  computation. 

What  is  the  size  of  the  matrix  A,  relative  to  the  memory  size  and  speed 
of  a  given  computer? 

If  we  are  solving  a  linear  equation  system  Ax  =  b,  do  we  have  many 
different  right-hand  sizes  b,  or  just  one?  Do  we  have  many  different  mat-rice 
that  are  close  together,  or  do  we  have  just  one  A?  Are  the  elements  of  A  ar 
precise  Mathematical  l  numbers  (for  example,  integers),  or  are  they  physical 
numbers  subject  to  uncertainty?  Any  uncertainty  in  A  and  b  leads  to 
uncertainty  in  the  definition  of  x  as  the  solution  of  Ax  =  b.  What  x  do-- 
the  problem's  proposer  want  to  see?  Even  when  A  and  b  are  mathematical 
numbers,  the  solution  x  is  normally  not  representable  as  a  finite-precision 
member  in  the  computer's  number  base.  Of  the  various  approximate  answers 
x  which  might  be  obtained,  what  is  the  proposer's  desire?  For  example,  does 
he  want  ||x  -  A~H>H  to  be  small,  where  A"H  is  the  true  answer?  Or  would  th 
proposer  settle  for  an  x  such  that  || Ax  -  b||  is  small?  For  each  case:  whir 
norm,  and  how  small? 

Most  proposers  of  linear  equation  systems  haven't  considered  these 
questions,  and  look  to  the  numerical  analyst  to  explain  the  possibilities  and 


select  the  options. 
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If  a  proposer  requests  the  inverse  matrix  A”1  ,  it  is  usually  worth  finding 
out  why.  Frequently  he  merely  wishes  a  convenient  way  to  solve  Ax  =  c  for 
an  arbitrary  vector  c.  Having  A  ^  stored  away,  the  proposer  expects  to 
obtain  the  solution  x  in  the  form  A~^c,  for  any  new  c  that  comes  along. 

It  should  be  pointed  out  that  there  are  other  ways  to  obtain  A”^c  for  new 
vectors  c,  ways  that  require  no  more  storage  and  take  no  longer  for  the  same 
accuracy,  than  the  multiplication  of  A-'*'  by  c.  Because  of  these  facts, 
the  computation  of  A  may  frequently  be  dispensed  with.  However,  certain 
statistical  applications  really  do  require  knowledge  of  at  least  the  diagonal 
elements  of  A  1. 

The  eigenvalue  problem  c  for  symmetric  matrices  A  can  require  finding 
all  the  eigenvalues,  or  only  a  few.  It  matters  a  good  deal  whether  or  not  the 
corresponding  eigenvectors  are  needed.  If  a  complete  set  of  eigenvectors  is 
needed,  is  it  important  that  they  be  orthogonal  to  each  other?  Getting  orthog¬ 
onal  eigenvectors  corresponding  to  multiple  eigenvalues  is  far  more  difficult 
than  just  getting  eigenvalues. 

In  the  eigenvalue  problem  d  for  nonsymmetric  matrices  A,  one  has 
similar  choices:  do  we  want  all  eigenvalues,  or  just  sane?  Do  we  want  column 
eigenvectors?  Do  we  want  row  eigenvectors?  Both?  But  then  comes  a  new  choice. 
If  some  eigenvalues  are  multiple  and  correspond  to  a  nonlinear  elementary  divisor, 
what  vectors  does  the  proposer  want  to  see?  In  monographs  on  algebra  one  learns 
about  chains  of  principal  vectors  that  with  the  eigenvector  form  a  basis  for  the 
null  space  N  of  (A  -  \  I)  ,  where  k  is  an  eigenvalue  of  multiplicity  k 
with  an  elementary  divisor  of  degree  k.  These  principal  vectors  are  associated 
with  the  Jordan  canonical  form  of  A.  It  is  my  impression  that  a  proposer  who 
has  a  good  background  in  algebra  will  want  to  see  a  set  of  principal  vectors 
(they  are  not  unique).  But  these  principal  vectors  are  extremely  hard  to  compute, 
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partly  because  they  are  discontinuous  functions  of  the  data.  It  is  likely  that 
an  orthogonal  basis  for  the  nullspace  N  would  be  a  more  useful  set  of  vectors 
The  matter  seems  to  be  poorly  understood  by  problem  proposers  and  numerical 
analysts. 

Matrices  with  actual  multiple  eigenvalues  are  very  rare,  and  a  small 
computational  perturbation  of  these  will  normally  destroy  the  equality  of 
eigenvalues.  One  might  therefore  assume  that  we  need  not  be  concerned  in  prac¬ 
tice  with  what  to  do  about  them.  But,  in  fact,  the  bad  behavior  of  nonlinear 
divisors  carries  over  in  practice  to  a  surprisingly  large  set  of  neighboring 
matrices.  These  neighboring  matrices  have  distinct  eigenvalues,  but  the  k 
column  eigenvectors  are  so  nearly  linearly  dependent  that  they  cannot  be 
separated  in  a  normal  computation.  So  also  here  one  faces  the  problem  of  what 
vectors  to  give  the  proposer. 

In  a  least  squares  problem,  say  a  search  for  x  to  minimize  f(x)  =  || Ax  ~bj;o, 
does  the  proposer  really  want  a  minimum  of  f(x),  or  does  he  merely  wish  an  x 
that  gives  a  value  of  f(x)  fairly  close  to  the  minimum?  In  a  curve-fitting 
problem,  for  example,  one  can  often  get  a  surprisingly  good  fit  by  a  polynomial 
with  coefficients  very  different  from  those  of  the  minimizing  polynomial. 

In  all  of  the  above  computational  problems,  it  is  important  to  ascertain 
which  of  the  following  types  of  answers  the  problem  proposer  is  looking  for: 

a)  a  surmised  answer,  with  no  estimates  of  its  correctness; 

b)  some  answer,  together  with  some  sort  of  probabilistic  assertions 
about  its  correctness; 

c)  some  answer,  together  with  mathematically  provable  bounds  for  its 

error. 

Normally  it  is  more  expensive  to  obtain  b)  them  a),  and  still  more 
expensive  to  obtain  c). 
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It  is  not  obvious  which  of  the  above  types  of  answer  the  problem  proposer 
will  want.  Frequently  a)  is  quite  satisfactory.  The  physical  scientist  and 
engineer  frequently  have  their  own  checks  on  the  validity  of  an  answer,  and 
may  neither  need  nor  wish  the  mathematician's  rigorous  bounds.  They  may  recog¬ 
nize,  for  example,  that  the  mathematical  model  is  such  a  rough  approximation  to 
reality  that  mathematical  bounds  would  only  be  ludicrous.  When  mathematicians 
enter  the  practical  world  of  engineering,  the  rules  by  which  mathematics  is 
played  frequently  have  little  relevance.  Numerical  analysts  frequently  have 
trouble  deciding  when  to  play  the  game  according  to  mathematician's  rules 
and  when  to  play  it  like  engineers.  It  is,  of  course,  extremely  pleasant  to 
encounter  those  occasional  exanples  where  mathematically  provable  bounds  can  be 
found  that  are  Just  as  accurate  and  cheap  as  surmised  answers.  One  should 
never  cease  looking  for  such  miracles,  because  they  do  occur!  One  has  been 
just  reported  at  this  SIAM  Symposium;  see  ?0x,  Henrici,  and  Moler  [26]. 

4.  Nature  of  computer  hardware  and  software.  The  character  of  achievable 
solutions  to  the  computational  problems  of  linear  algebra  is  greatly  influenced 
by  the  nature  of  the  computing  systems  available  to  us.  It  is  customary  to 
separate  computer  systems  along  the  following  lines: 

a)  Computer  hardware--the  nature  of  the  electronic  circuitry  of  a 
computer; 

b)  Computer  languages — the  languages  in  which  are  described  algorithms 
for  the  solution  of  a  given  problem  on  a  given  computer; 

c)  Computer  software — the  programs  which  make  it  possible  for  a  computer 
actually  to  perform  the  algorithms  described  in  the  computer  language. 

In  looking  at  computer  hardware  for  computations  in  linear  algebra  one 
wants  to  know  what  precision  is  available  for  computation--how  many  digits  are 
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in  the  significant!  of  the  floating-point  operands,  and  to  what  base?  One  :s 
also  interested  in  the  cost  and  speed  of  double-precision  operations.  In 
matrix  algebra  work  the  critical  operation  is  frequently  the  computation  of 
a  rounded  single-precision  approximation  to  the  double- precision  inner  product, 
of  two  vectors  whose  components  are  single-precision  floating-point  numbers. 

The  speed  and  cost  of  this  inner  product  are* very  important. 

One  wonders  whether  the  hardware  rounds  the  result  of  an  arithmetic 
operations,  or  whether  it  is  chopped  off.  Best  of  all  is  a  system  that  lets 
the  programmer  decide  when  to  round  and  when  to  chop. 

What  happens  when  the  result  of  an  arithmetic  operation  exceeds  the 
capacity  of  the  floating-point  system?  Are  there  "traps"  which  make  it  possible  for 
the  system  to  detect  overflow  or  underflow?  Can  these  traps  be  by-passed, 

(turned  off)  by  the  programmer?  When  an  overflow  or  underflow  is  detected,  is 
all  essential  information  recoverable,  so  that  the  solution  can  continue?  Or 
are  vital  bits  of  information  irretrievably  lost? 

What  is  the  exact  nature  of  the  arithmetic  operations  in  the  machine?  IT 
one  is  to  prove  theorems  about  the  behavior  of  a  computation,  one  needs  certain 
properties  of  the  arithmetic..  Because  of  the  rounding  of  the  machine,  it  is 
well  known  that  addition  and  multiplication  are  not  associative,  nor  are  they 
distributive.  Nevertheless,  one  can  do  surprisingly  good  analysis,  provided 
only  that  the  arithmetic  is  monotonic . 

By  multiplication  being  monotonic,  we  mean,  for  example,  that  if 
0  <  a  <  b  and  0  <  c.  then  a  X  c  <  b  X  c  .  Such  properties  seem  elementa.. . 
out  they  are  extremely  helpful.  And  they  are  surprisingly  often  absent'. 

It  must  be  noted  that  apparently  minor  changes  in  the  hardware  of  the 
arithmetic  circuitry  can  make  surprisingly  large  differences  in  the  behavior 
of  the  algorithms. 
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A  great  many  computer  languages  have  been  devised  for  the  description  of 
scientific  algorithms.  These  range  from  the  very  elementary  codes  for  Turing 
machines,  through  the  machine  codes  of  computers,  to  various  algebraic  languages 
like  the  forms  of  Fortran,  Algol,  and  PL/l.  All  these  languages  are  equiva¬ 
lent,  in  the  sense  that  the  class  of  representable  algorithms  ia  the  same  for 
all  of  them.  The  languages  differ  only  in  regard  to  human  convenience  and  in 
the  compilation  problems  they  create.  Can  one  conveniently  represent  such  a 
data  structure  as  a  triangulr  ■  matrix  in  a  certain  language?  In  typical 
languages  like  Algol  or  Fortran,  one  must  choose  between  representing  it  as 
part  of  a  much  larger  square  matrix,  on  the  one  hand,  or  as  an  artificially 
created  one -dimensional  array,  on  the  other.  The  former  choice  is  humanly 
convenient  and  wastes  space;  the  latter  choice  saves  the  computer  time  and 
space,  at  the  cost  of  confusing  the  human. 

Most  matrix  algorithms  have  "inner  loops"  where  most  of  the  computing 
time  is  spent.  If  only  this  inner  loop  is  programmed  very  efficiently  in 
machine  code,  the  program  will  run  very  rapidly.  It  scarcely  matters  how  the 
rest  of  the  algorithm  is  programmed.  Hence  a  very  important  question  for  any 
algebraic  language  is  whether  it  .1$  -■  easy  to  incorporate  pieces  of  machine  code 
into  them.  Perhaps  the  question  is  more  appropriately  addressed  to  the  software 
system  that  translates  the  algebraic  language  into  machine  code. 

Another  important  property  of  a  computer  language  is  its  readability  by 
human  beings.  If  the  algorithm  is  correctly  written,  a  computer  will  (practically 
always  read  it  correctly.  But  the  practical  use  of  the  algorithm  depends  on 
the  ability  of  human  beings  to  comprehend  it,  adapt  it  to  other  uses,  improve  it 
in  the  light  of  recent  discoveries,  and  so  on.  The  human  readability  of  existing 
languages  differs  a  great  deal. 
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The  most  important  software  programs  for  the  scientific  computer  user 
are  the  monitors  and  the  compilers.  The  compilers  are  vast  symbol-manipulation 
programs  that  translate  an  algorithm  from,  say,  Fortran  to  the  machine  code  of 
a  given  computer.  Compilers  should  be  distinguished  from  the  languages  they 
translate,  and  yet  of  course  compilers  and  languages  influence  each  other. 
Compilers  differ  greatly  in  speed,  in  the  optimality  of  the  machine  code 
produced  in  the  translation,  and  in  the  diagnostic  facilities  offered. 

As  we  noted  above,  it  is  important  that  compilers  be  able  to  accept 
pieces  of  algorithms  written  in  machine  code,  and  incorporate  them  into  a  program 
otherwise  written  in  an  algebraic  language.  For  matrix  work,  the  ability  to 
compile  fast  codes  for  iterative  loops  (the  for  statment  of  Algol)  is  very 
important . 

Most  compilers  are  now  imbedded  in  control  programs  variously  called 
master  control  programs,  monitor  systems,  or  operating  systems.  These  monitor 
systems  generally  retain  ultimate  control  of  a  computer,  preventing  a  possibly 
erroneous  user  program  from  consuming  vast  amount  of  unwanted  time,  or  from 
damaging  the  monitor  system  or  other  persons'  programs  by  illegal  assignments 
Also,  the  monitor  systems  generally  recover  control  of  the  machine  in  case  of 
overflow  or  underflow.  This  is  a  point  of  much  interest  to  writers  of  linear 
algebra  programs.  In  case  of  overflow  or  underflow,  what  happens  next?  Can 
the  linear  algebra  program  recover  control  of  the  computer  and  repair  the  dama 
done  by  the  overflow  or  underflow?  (This  assumes  that  the  hardware  retains 
necessary  information. )  Or  does  the  monitor  system  take  over  the  machine  and 
ruthlessly  flush  the  offending  program  from  the  machine?  If  the  latter  occurs, 
then  extra  time  must  be  taken  in  each  program  to  make  sure  that  overflow  or 
underflow  cannot  occur. 


+~nm*  -Tt  - 
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5.  The  state  of  the  art,  1953  and  now. 


It  is  safe  to  say  that  matrix 


computation  has  passed  well  beyond  the  stage  where  an  amateur  is  likely  to 
think  of  computing  methods  which  can  compete  with  the  better  known  methods. 
Certainly  one  cannot  learn  theoretical  linear  algebra  ?>nd  an  algebraic 
programming  language,  and  nothing  else,  and  start  writing  programs  which  will 
perform  acceptably  by  today's  standards.  There  is  simply  too  much  hard-earned 
experience  behind  the  better  algorithms,  and  yet  this  experience  is  hardly 
mentioned  in  mathematical  textbooks  of  linear  algebra. 

The  amount  of  literature  on  matrix  computations  is  staggering.  In  620 
pages,  Faddeev  and  Faddeeva  [8]  record  a  pretty  complete  account  of  computational 
methods  up  to  around  1958.  In  662  pages,  Wilkinson  [61]  gives  most  of  what 
is  known  about  computing  eigenvalues  of  dense,  stored  matrices  (both  symmetric 
and  unsymmetric),  with  error  bounds  for  many  algorithms.  There  is  very  little 
overlap  between  the  two  books,  because  Wilkinson  and  a  few  contemporaries 
created  most  of  the  material  in  his  book  in  the  years  after  1958.  No  one  could 
possibly  start  research  in  the  numerical  mathematics  of  linear  algebra  without  a 
thorough  knowledge  of  the  relevant  material  in  these  books. 

It  is  perhaps  instructive  to  examine  the  state  of  matrix  computation 
in  1955,  when  the  author  wrote  a  survey  [10]  of  methods  for  solving  linear 
systems  at  the  Institute  for  Numerical  Analysis  of  the  National  Bureau  of 
Standards,  Los  Angeles.  We  were  amateurs.  For  dense,  stored  matrices  we 
knew  Gaussian  elimination,  of  course.  We  knew  that  it  sometimes  produced 
ouite  poor  results.  We  weren't  always  sure  why.  We  debated  endlessly  about 
how  to  pick  pivots  for  the  elimination,  without  settling  it.  The  debate  still 
continues, but  now  mainly  among  persons  who  don't  understand  that  the  main  lines 
of  the  answer  have  been  settled.  Because  of  misunderstood  difficulties  with 
Gaussian  elimination,  we  searched  for  other  methods  which  might  do  better. 
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The  conjugate-gradient  method  had  been  devised  for  sparse  matrices  by  Lanczos 
(36],  and  Hestenes  and  Stiefel  (27].  In  [lO]  I  guessed  that  it  might  also  prevail 
for  dense,  stored  matrices,  despite  the  extra  time  it  would  require,  because  we 
understood  how  to  use  higher  precision  to  make  the  conjugate-gradient  method  work 
well.  We  did  not  realize  that  the  same  higher  precision  and  a  proper  pivotal 
strategy  would  make  Gaussian  elimination  work.  We  were  not  quite  aware  of  the 
extent  of  problems  of  ill  conditioning  of  matrices. 

The  only  analysis  available  to  us  was  the  monumental  work  of  von  Neumann 

and  Goldstine  [4l,  20].  They  avoided  the  pivoting  problem  by  reducing  any 

regular  linear  equation  system  Ax  =  b  to  the  positive  definite  system 
T  T 

A  A  x  =  A  b.  We  knew  that  this  normalization  of  the  problem  was  costly  in  time 
and  worsened  the  condition  of  the  problem.  Von  Neumann  and  Goldstine  presented 
guaranteed  error  bounds  for  the  solution;  actually  observed  errors  were  found 
to  be  perhaps  100  times  smaller  in  reasonable  cases.  The  form  of  the  error 
analysis  was  a  direct  comparison  of  machine  arithmetic  with  exact  operations. 

The  nonassociativity  and  nondistributivity  of  machine  arithmetic  made  the 
analysis  extremely  difficult.  In  any  case,  it  could  only  handle  scaled  fixed- 
point  arithmetic.  Because  of  the  size  of  their  error  bounds,  von  Neumann  and 
Goldstine  were  unnecessarily  pessimistic  about  the  possibility  of  inverting 
general  matrices  of  orders  over  15  on  machines  with  the  27-bit  precision  of  the 
IBM  7090  series. 

For  the  eigenvalue  problems,  things  were  In  much  worse  state.  We  had  the 
power  method  with  matrix  deflation.  While  reasonably*  satisfactory  for  a  few 
dominant  roots,  its  general  application  requires  intuition  and  luck,  and  defies 
a  complete  algorithmization.  For  dense,  stored  symmetric  matrices  we  had  the 
1846  method  of  Jacobi  (3l]>  rediscovered  and  analyzed  by  Goldstine,  Murray 
and  von  Neumann  [19L  and  it  was  quite  satisi'actory.  Givens  was  writing  up 
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his  newly  discovered  method,  maybe  7  to  9  times  faster  than  Jacobi's  and 
a  basic  step  toward  currently  used  methods. 

For  nonsymmetric  matrices,  things  were  ghastly.  If  the  power  method 
wouldn't  work,  we  had  practically  no  alternatives.  We  could  search  for  zeros 
of  det(A  -  zl)  in  some  manner  or  another.  We  bravely  tried  methods  for  deter¬ 
mining  the  characteristic  polynomial,  as  described  in  Faddeeva  [9]*  and  found 
them  to  be  hopeless.  It  was  almost  unbelievable,  how  badly  the  standard 
methods  for  n  =  4  would  perform  for  n  =  10.  Lanczos  was  advocating  his 
new  method  of  finite  iterations,  which  became  the  source  of  modern  methods 
in  a  later  line  of  development  through  the  Stiefel  and  Rutishauser  QD-algorithm, 
(see  Rutishauser  [50]  and  Henrici  [ 2 5 ]) >  the  LR-algorithm  of  Rutishauser  [5l]> 
and  the  QR  algorithm  of  Francis  [l5,  l6]  and  Kublanovskaja  [35].  However,  the 
original  Lanczos  method  needed  careful  management,  because  the  raw  results 
were  often  poor. 

6.  The  linear  equations  problem.  For  large,  sparse  matrices,  like  those 
arising  in  finite-difference  approximations  to  partial  differential  equations, 
there  is  a  whole  special  literature.  See  Varga  [57L  Forsythe  and  Wasow  [13], 
the  work  of  David  Young,  Jim  Douglas  Jr.,  Stiefel,  and  many  others.  The  methods 
seem  to  depend  for  their  success  on  the  nature  of  the  continuous  problem  being 
approximated  Because  the  matrices  are  sparse,  the  prevailing  methods  are 
iterative.  I  shall  omit  further  discussion  of  them,  and  confine  attention  to 
dense,  stored  matrices. 

For  a  general  matrix  A,  the  solution  of  the  linear  system  Ax  =  b  by 
Gaussian  elimination  requires  r?h  +  0(n2)  multiplications,  and  the  same 
number  of  additions.  Recently  Klyuyev  and  Kokovkin-Shcherbak  [3^]  proved  that 
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no  method  using  ration'll  operations  for  general  A,  b  can  take  fewer  operations. 
This  result  had  long  been  believed  but  not  proved.  The  result  has  two 
consequences: 

(i)  Gaussian  elimination  is  likely  to  remain  the  method  of  enoice  for 
solving  dense  linear  systems,  when  it  works,  because  it  is  as  fast  as  any. 

(ii)  The  solution  of  a  linear  system  of  large  order  n  is  going  to 
require  a  very  substantial  amount  of  computing  time,  cat  least  for  serial 
computers.  For  n  =  1000,  we  have  l/3  x  10^  multiplications  and  additions. 

If  we  can  multiply  and  add  in  10  microseconds,  we  need  3333  seconds,  or  about 
an  hour  of  computation.  In  fact,  there  is  some  overhead  also,  and  on  an  IBM 
7094  (Model  II )  the  solution  would  take  over  2  hours.  However,  the  storage  cf 
the  million  elements  of  data  requires  extensive  use  of  some  bulk  storage  like 
tapes  or  disks,  as  only  some  20,000  elements  or  so  can  be  kept  in  the  current 
32,000-word  core  storage.  The  very  numerous  transfers  of  matrix  elements  from 
core  to  magnetic  tapes  appear  likely  to  wear  out  the  tapes  before  the  solution 
can  be  obtained,  according  to  certain  tests  made  at  Stanford;  I  know  of  no 
comparable  experience  with  magnetic  disks  or  other  form  of  bulk  storage. 

As  a  result,  we  cannot  consider  order  n  =  1009  to  represent  a  practical 
linear  equations  problem,  but  we  will  undoubtedly  soon  be  able  to  do  i+.  regularly 
for  perhaps  $500. 

The  case  n  =  100  is  now  easy  and  costs  around  $1  on  an  IBM  709*+.  The 
case  n  =•  10,000  is  likely  not  to  be  accessible  for  a  long  time,  and  it.  would 
take  over  2000  hours  now  on  an  IBM  7094. 

There  is  beginning  to  be  serious  consideration  of  computers  with  a 
substantial  amount  of  parallel  operation,  so  that  perhaps  much  of  the  solution 
of  n  linear  system  could  be  done  simultaneously.  Preliminary  studies  make  it 
clear  that  the  solution  of  a  linear  system  could  very  easily  make  use  of  para J  h  1 
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computation,  if  it  should  prove  worth  while.  Apparently  only  C(n)  operation 
times  would  be  needed  for  solving  a  linear  system,  if  one  had  a  sufficiently 
large  amount  of  parallel  arithmetic  circuitry. 

7.  Inherent  inaccuracy  in  solutions  of  linear  systems.  Given  a  non¬ 
singular  matrix  A  and  a  nonzero  b,  let  x  be  the  solution  of  Ax  -  b. 

Suppose  A  and  b  are  subject  to  uncertainty.  What  is  the  resultant  uncer¬ 
tainty  in  the  solution  x? 

For  any  column  vector  y  of  order  n,  define  ||y||  to  be  the 
euclidean  length  of  y: 

llyll  =  llyll2  =  J  +  A.  +  ’  ‘  *  +  yn  • 

For  any  n-by-n  square  matrix  A,  define  the  spectral  norm  ||A||  by 

||a||  =  max  ||Ax||  . 

Ilxll-l 

These  functions  ||...||  give  useful  measures  of  the  size  of  vectors  and  matrices, 
respectively. 

For  a  nonsingular  matrix  A,  define  the  condition  of  A,  cond(A) 
by  the  relation 

cond(A)  =  I!a||  -  ||A_1|| 

The  concept  of  condition  of  a  matrix  seems  to  have  been  introduced  by  Turing  [55] » 
and  studied  extensively  by  Todd  ([55l  and  some  later  papers)  and  many  others. 

One  of  the  main  uses  of  the  concept  of  condition  lies  in  answering  the 
question  posed  at  the  start  of  this  section.  Suppose  that  A  is  known  exactly, 
but  that  b  is  subject  to  uncertainty.  Let  x  +  6X  solve  the  system  with  matrix 
A  and  right-hand  side  b  +  db.  Then 
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(3) 


A(x  +  &x)  =  b  +  db  ; 

x  +  6X  =  A’H  +  A_1db  ; 
&x  =  A~  ^db  ; 

IM  <  ||A_1i|  •  |M|  . 

Since  A  x  =  b,  we  have 


(4) 


Ml  =  li Ax |  <  ||a||  •  ||x| 


Dividing  (3)  by  (4),  we  have 


UMll  <  ||A||  .  HA-1!!  JIM 

11*11  "  !lt!l 


or 

(5>  M  <  Cond(A)lM  . 

11*11  "  llb'l 

Inequality  (5)  shows  that  the  relative  uncertainty  in  x  is  bounded  by 
eond(A)  times  the  relative  uncertainty  in  b.  The  bound  in  (5)  is  attainable, 
for  any  nonsingular  A  and  nonzero  b.  This  is  easy  to  see,  if  we  perform  a 
change  of  coordinates  in  which  A  takes  a  diagonal  form. 

An  a  .linear  transformation,  A  takes  vectors  x  into  vectors  b. 

A  fundamentally  important,  but  too  little  known  theorem  states  that  by  a  certain 
orthogonal  change  of  coordinates  in  the  space  of  x,  and  by  another  orthogonal 
change  of  coordinates  in  the  space  of  b,  the  matrix  A  can  be  put  in  the  diagonal 
form 
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Here  the  positive  numbers  >  \l^  >  . . .  >  M-n  are  called  the  singular  values 
of  A.  Moreover, 


Finally,  the  orthogonal  transformations  do  not  change  the  norms  of  x  and  b. 


We  have 
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ildbll 

Ml 


cond(A)  IM 

IN 


For  these  vectors, 

IM 

tun  % 

The  last  line  shows  that  (5)  is  an  equality  in  this  case,  as  we  promised 
to  prove. 

Although  (5)  is  only  an  exact  equality  under  exceptional  conditions, 
it  is  usually  rather  close  to  equality,  and  in  the  following  we  assume  approxi¬ 
mate  equality. 

If  cond(A)  =  IQ1*  ,  and  if  b  is  known  to  be  correct  only  to  10  decimals, 
then  x  can  be  known  only  to  10  -  p  decimals.  Now  p  can  range  anywhere  from 
0  to  oo  .  The  only  hope  of  having  any  significance  to  x  in  a  10-decimal  computing 
system  is  that,  roughly, 

cond(A)  •  10" 10  <  |  . 

In  a  base -3  computer  with  t  significant  digits,  we  ro  .ghly  need 

cond(A)  ^ 

in  order  to  have  any  significance  to  a  solution. 

Remember  that  all  statements  in  this  section  are  independent  of  any 
method  of  solving  a  system  Ax  =  b.  They  are  statements  about  errors  in  x 
which  are  inherent  in  the  uncertainty  in  the  data. 

If  A  is  subject  to  a  change  dA,  and  b  is  known  exactly,  then 
an  inequality  analogous  to  (5)  is  the  following: 

■  JIM  <  cond(A)  *  M- 

||x  +  bx||  ||a|| 

If  ;jcxij  is  small,  compared  with  ||x||,  then  we  may  safely  consider  the  left-hand 
si  e  of  (<>)  as  a  relative  error  in  x. 
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8.  Accuracy  achievable  with  Gaussian  elimination.  I  assume  that  the 
reader  knows  what  Gaussian  elimination  is,  as  a  method  of  solving  linear  eaua- 
tion  systems.  The  main  strategic  decision  facing  the  designer  of  the  algorithm 
is  the  choice  of  a  unique  pivot  element  for  each  of  the  n-1  stages  in  which 
a  variable  is  eliminated  from  the  remaining  equations.  There  are  two  main 
strategies  discussed: 

(i)  complete  pivoting,  in  which  at  each  stage  one  selects  as  a  pivot 

some  element  a.  .  of  maximum  absolute  value  among  all  the  remaining  elements 
1  j  0 

of  the  matrix. 

(ii)  partial  pivoting,  in  which  at  each  stage  one  selects  as  a  pivot  some- 

element  a.  of  maximum  absolute  value  among  the  first  column  of  the  remainin’: 

1  >  J 

elements  of  the  matrix. 

Thus,  in  the  first  stage  complete  pivoting  would  search  the  whole  matrix 
A  for  an  element  maximal  in  absolute  value,  whereas  partial  pivoting  would 
search  only  the  first  column. 

Some  special  classes  of  matrices  permit  elimination  to  proceed 
successfully  without  any  search  for  pivoting — for  example,  positive  definite 
symmetric  matrices.  But  generally,  pivotal  searching  is  essential  to  guarantee 
success.  The  following’  simple  example  illustrates  the  disaster  possible  in  not 
searching  for  a  pivot.  Consider  a  3-digit  floating-decimal  machine. 

The  system  is 


{ 


.0001  x  +  1.00  y  =  1.00 

1.00  x  +  1.00  y  =  2.00  . 


The  true  solution,  rounded  to  five  decimals,  is  x  =  1.00010,  y  =  .99990. 

If  one  accepts  the  element  .0001  as  a  pivot,  the  elimination  of  x 

from  the  second  equation  yields  the  equation 

-  10000  y  =  -10000. 

Backsolving,  we  find  that  y  -  1.00,  whence  x  =  0.00,  a  clear  disaster. 
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On  the  other  hand,  partial  pivoting  would  select  the  element 
,  =  1.00  as  the  pivot.  Elimination  of  x  from  the  first  equation  yields 

> 

the  equation 

1.00  y  =  1.00  . 

Backsolving,  we  get  y  =  1.00  and  then  x  =  1.00,  with  obvious  success. 

We  shall  now  assume  that  we  are  dealing  with  a  t-digit  base-2  floating¬ 
point  computer.  Rather  than  discuss  the  solution  of  a  linear  system,  we  shall 
consider  the  computation  of  the  inverse  A*1  of  a  given  matrix.  We  wish  to 
state  the  rounding  error  bounds  that' have  been  proved  for  Gaussian  elimination. 

Wilkinson  [58]  assumes  a  complete  pivotal  strategy,  and  that  the  matrix  A 
is  reasonably  scaled  at  the  start  and  at  all  intermediate  stages  (see  Sec.  10 

ui 

yields  a  matrix  X  such  that 


for  more  about  scaling).  Then,  if  all  ; 


|  <  1,  a  certain  Gaussian  algorithm 


(7)  <  (0.8)2-*  n7/2g(n)  Ia'1!  . 

Here  g(n)  is  the  maximum  of  all  elements  of  the  successive  matrices  found 
during  the  elimination. 

To  express  the  result  (7)  in  a  form  to  be  compared  with  those  of  Sec.  7» 

•  l/2 

we  note  that  1  <  |jA||  <  n,  so  that  we  expect  that  |l A}[  =  n  '  .  Then  we  have 

roughly 

(8)  iX  ~.A~1J1  <  n3  •  2"*  cond(A)  g(n)  . 

Ia-4  - 

What  kind  of  bound  can  we  give  for  g(n)?  This  turns  out  to  be  an  open 
question.  The  best  known  result  is  approximately 

(9)  g(n)  <  1.8  n^AUog  n. 


22 


On  the  other  hand,  for  I  j  real  mat  r  i  r •  t,  r»ver  xam'ir.ed  it  has  always  been 
observed  that 

g(n)  <  n  . 

The  last  bound  is  attained  for  unboundedly  large  n  by  matrices  related  to  the 
H  idamard  matrices.  For  most  matrices  one  even  observes  that 

(10)  g(n)  <  8  . 

Tornheim  [5*+J  has  found  complex  matrices  A  of  unboundedly  large  n  for 

if  .ich 

g(n)  =  5.1  n  . 

It  would  be  most  desirable  to  have  a  good  bound  for  g(n),  so  that  (7) 
could  be  turned  into  a  good  a  priori  error  bound  for  the  computation  of  A  1  . 
Naturally,  for  any  particular  matrix  A,  g(n)  is  easily  observed  in 

the  course  of  the  elimination,  so  that  in  any  event  (7)  becomes  an  a  posteriori 

« 

error  bound.  However,  still  better  error  bounds  can  be  given  a  posteriori,  as 
will  be  shown  in  Sec.  9. 

Wilkinson's  proof  of  (7)  in  [58]  is  reasonably  short.  It  makes  use  of 
'averse  rounding  error  analysis,  which  we  shall  mention  again  in  Sec.  11.  It  is 
'..istructi ve  to  compare  (8)  with  (6),  even  though  one  deals  with  inverses  and  one 
..'.th  linear  systems.  The  factor  2“^  is  essentially  the  inherent  uncertainty 
’evel  of  the  data,  and  should  be  equated  to  ||dA||/|| A|| ..  Then  the  bound  in  (8) 
is  larger  than  that  in  (6)  by  the  factor  n  g(n),  Taking  into  account  the  empir¬ 
ical  result  (lO)  that  g(n)  <  8  for  most  real  matrices,  we  then  interpret  (8) 

•  3  saying  that  the  computed  matrix  X  generally  differs  from  the  true  inverse 
\  in  relative  terms  by  no  more  than  8n  times  the  error  inherent  in  the 
problem.  Thus  simple  Gaussian  elimination  is  reasonably  good  at  keeping  the 
rounding  error  bound  under  control,  for  modest  values  of  n.  Much  better  results 
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can  be  achieved  with  some  devices  to  be  mentioned  in  Sec.  9* 

The  bound  corresponding  to  (7)  given  by  von  Neumann  and  Goldstine  [^l] 
was 


(11) 


l|A-l|| 


<  (5.3  +  1 4.6  l|A||2)  2-t  n2  ||A-1||2 


15  n2  2“t[cond(A)]2 


2  T  T 

The  factor  [cond(A)]  arose  from  solving  A  A  x  =  A  b,  rather  than  Ax  -  b. 
The  proof  of  (ll)  was  an  order  of  magnitude  more  difficult  and  tedious  than 
the  proof  of  (7). 


9.  More  accurate  solutions.  Suppose  that  A  is  given  as  single¬ 
precision  data,  and  that  we  wish  to  get  solutions  guaranteed  to  be  more  accurate 
than  the  above  bounds  would  indicate.  How  shall  we  proceed?  The  most  obvious 

choice  is  to  perform  all  calculations  in  double-precision.  Roughly  speaking, 

-2t. 

then  t  is  replaced  by  2t  in  the  above  error  bounds,  and,  since  2  is 
so  very  much  smaller  than  2_t,  we  gain  many  orders  of  magnitude  in  acci.racy. 

The  cost  in  computing  time  varies  among  different  machines,  but  is  only  a  factor 
of  four  on  the  IBM  709^.  The  cost  in  storage  is  greater,  since  we  must  double 
the  storage  reserved  for  the  developing  matrix. 

Where  the  time  and  storage  costs  are  too  high  to  justify  complete  double 
precision,  it  is  possible  to  make  a  very  substantial  gain  by  a  much  ir.01  e  limited 
use  of  double  precision.  Most  of  the  operations  in  Gaussian  elimination  can 
be  phrased  as  inner  products  of  vectors  of  single-precision  numbers.  Gn  many 
machines  it  is  possible  to  accumulate  such  an  inner  product  in  double  precision, 

'  nd  then  round  it  off  to  single  precision  before  storing  away  the  result. 

The  result  of  this  accumulation  is  to  reduce  the  maximum  rounding  error  of  an 

2k 


inner  product  by  a  factor  of  n.  The  total  effect  turns  out  to  be  to  reduce  the 

5/2 

round-off  error  bound  by  a  factor  of  n  '  .  Thus,  instead  of  the  result  (7)> 
an  elimination  with  pivoting  and  accumulation  produces  an  approximate  inverse  X 
such  that 

(12)  N*  -77—^  <  3.3n  2'1  HA-1!!  g(n)  , 

under  certain  additional  hypotheses.  See  Wilkinson  [6l,  p.  253].  The  gain 
5/2 

of  the  factor  n  is  very  substantial,  although  experience  shows  that .  the 
actual  errors  in  single-precision  computation  are  usually  rather  less  than  the  bounds. 

One  theoretical  disadvantage  of  the  complete  pivoting  strategy  is  that 
it  does  not  mix  well  with  the  accumulation  of  inner  products.  When  products 
are  accumulated,  one  almost  always  uses  a  partial  pivotal  strategy,  and  accepts 
the  theoretical  possibility  that  pivots  can  grow  very  large. 

A  third  and  the  most  successful  approach  to  increasing  the  accuracy  of 
solutions  of  dense,  stored  linear  systems  is  the  so-called  method  of  iterative 
improvement .  By  this  method,  if  the  matrix  A  is  not  too  ill-conditioned,  one 
in  practice  gets  solutions  which  are  the  correctly  rounded  approximations  to  the 
truj  answers.  We  will  now  describe  this  development. 

Suppose  that  by  Gaussian  elimination  one  has  achieved  a  first  approximate 
solution  Xq  of  the  linear  system  Ax  =  b.  The  next  step  is  to  form  the  residual 
vector  rQ  =  b  -AxQ.  If  xQ  were  the  exact  solution  of  the  system,  we  would 
have  rQ  =  0  ,  the  null  vector.  If  not,  we  solve  the  new  linear  system  Ay  =  rQ  , 
to  obtain  a  vector  y^.  Let  =  xQ  +  y^. 

The  process  is  repeated  iteratively.  I.e.,  for  k  =  0,  1,  2,  ... 

,  solve  the  system  Ay  =  r^  to  obtain  a 
vector  yR+1,  and  then  form  xk+1  =  \  +  Yk+1  • 


we  form  the  residual  r, 


-b  -  *** 
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Under  suitable  hypotheses  to  be  specified  below,  the  sequence  x 

K. 

-l 

converges  to  the  true  solution  A  b  of  the  system  Ax  -  b. 

Several  matters  need  to  be  clarified  in  this  algorithm.  First,  it  appears 
to  involve  a  great  deal  of  work  to  solve  systems  of  the  form  Ay  r^  for  many 
values  of  k.  In  fact,  this  is  not  so.  Gaussian  elimination  to  solve  a 
system  Ax  =  b  involves  three  distinguishable  stages: 

(i)  Triangularization  of  the  matrix  A  by  elementary  row  : rar.sfcrma*  :ens 

(ii)  Application  of  the  same  row  transformations  to  the  right-hand  sid'-'  to 

(iii)  Solution  of  the  triangular  system  by  back-substitution. 

V 

It  turns  out  that  stage  (i)  requires  approximately  n/5  multiplies!  ions 

and  additions,  but  that  stages  (ii)  and  (iii)  together  requre  ^nly  apprcxima*  ely 
2 

n  multiplications  and  additions.  Stage  (i)  need  be  done  only  once  for  axl 

the  systems  Ay  -  r^.  If  the  multipliers  defining  the  row  transformations  are 

saved,  stages  (11)  and  (iii)  can  be  done  rapidlj  for  each  new  system  Ay  r^ 

m  turn.  As  o  result,  it  is  found  that  b  sufficiently  long  sequence  of  vectors 

x  con  usually  be  computed  in  something  like  only  20  per  cent  more  time  1  nan 

the  computation  of  the  first  solution  x^. 

It  is  absolutely  essential  that  each  residual  vector  r  oe  comp^’ed  •  c 
_ K _ 

nigh  precision,  '"his  is  normally  done  by  a  do-ble-preci sior.  ^ccum-ia'  ion  of 
inner  prod-cts,  followed  by  rounding  of  the  answer  to  single-precision  lloa* ing- 
pom*,  form.  If  r^  is  computed  with  merely  a  singxe-precisicn  inner  pr  d.ct  . 
it  will  have  rounding  errors  of  several  units  m  ’he  least  significant  e.g.’s 
of  x,  '.’.'hen  the  inequality  (5),  which  :s  an  approximate  eq-niity  in  practice, 


tells  us  that, 

significan4  digit 

accuracy  m  x. 

k 


will  oe  wrong  by  several  times  cond(A)  in  its  leas; 

/x  US 

.  Since  cond(AJ  may  well  be  10  or  10  ,  the  result  an' 
is  very  low  and,  in  fact,  is  itself  almost  as  acc..r-te 


is  '<ny  succeeding  x  . 

K. 
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*'hf  following  the  ore  k.  gives  4  n*  ens  i  s  .f  the  above  me-,  hod  of 
iterative  improvement: 

Tneorem  Let  the  ratrix  A  nave  ‘  he  property  that 

(15)  (0  8)  2"'  n7/  '  c(r)  Ha*1!!  <  *  , 


_et  the  above  algorithm  be  carried  ou-.,  v:Th  each  system 

_n  single-precision  base-2  float  ing-pcint  antnnietic,  but 

r  b  -  Ax,  and  x,  , .  -  x,  +  y  .  carried  o.t  without 
k  k  -  k+1  k  k+1  - 

'nen  , 

j|x  -  A  b||  -»  0,  as  k  -»  oi  , 

K 


Ay  -  r  being  solved 
K 

with  computations  of 
rOw.nd.tng  error. 


J.f  the  solution  cf  the  systerrs  Ay  -  r  were  done  with  acc -dilations  of 

K 

.  r.ner-products  in  double  precision,  tnen  the  left-hand  side  of  (l3)  could  be 

* ^placed  by  the  right-hand  side  of  (i 2) 

i.n  practice,  of  course,  r  ,s  comp^f-'d  oy  a  double-precision 

K. 

cumulation  of  inner  products,  and  x  ^  is  comp. ted  as  the  float mg-po.nt 

■  of  x.  and  y  As  a  resui4  ,  the  sequence  x  does  not  converge  tc 

K  Kr  1  K. 

b  m  the  mathematical  sense.  Instead,  x  is  observed  to  oecon.e  constant 

k 

a  value  wnich  is  normally  the  coirectly  rounded  single-precision  approximation 
A_1b 

In  the  actual  use  of  iterative  improvement,  one  does  not  -Sw.aliy  know 
advance  whether  or  not  hypothesis  f  13)  is  satisfied,  and  it.  cannot  be  con- 
.ied  afterwards  either.  Normal  practice  is  there  tore  to  rely  on  the  following 

■  -  .r: stic  result " 

Almost -theorem.  Let  the  above  algoritnm  be  carried  Ow.t ,  with  each 
.-..’’em  Av  r,  being  solved  by  4  he  same  version  of  Ga-ssian  elimination,  with 

— —  —  K.  * 

u  being  computed  by  a  double-precision  gccuivla*  ion  of  inner  products, 

-to  with  xk  +  ,  being  computed  as  the  float  ing- point  s-m  of  x^  and 
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If  for  k  >  all  vectors  x^  are  equal  to  some  single-precision 
vector  x  ,  then  x  is  the  correctly  rounded  single-precision  approximation 
to  A-1b  . 

This  almost-theorem  cannot  be  proved,  and,  indeed,  Kahan  [52]  has  an 
extremely  ingenious  counter-example.  However,  roost  computers  would  bet  their 
life  on  the  applicability  of  the  above  almost-theorem  in  any  practical  example, 
unless  Kahan  were  furnishing  the  problem! 

Normally,  when  cond(A)  gets  near  2^,  the  vectors  x^  obviously 
diverge.  Then  there  is  no  cure  except  to  increase  the  precision  with  which 
the  elimination  is  carried  out,  unless  scaling  A  will  help. 

The  usual  value  of  k_  is  3  or  4. 

u 


10.  Scaling  of  matrices.  One  matter  that  was  glossed  over  in  Sec.  8 
was  the  scaling  of  the  matrix  A  before  solving  a  system  Ax  =  b.  Alternate 
terms  for  scaling  are  preconditioning  and  equilibration.  Suppose  that  the 
2-by-2  numerical  example  of  Sec.  8  were  altered  by  multiplying  the  first 

100000 

2.00 

The  effect  of  the  scaling  is  to  make  10.0  the  larger  pivot  in  the  first  column. 
Then  elimination  of  x  from  the  second  equation  of  (l4)  in  5-digit  floating- 
decimal  arithmetic  will  result  in  a  new  second  equation 

-10000  y  =  -  10000  . 

Back  solution  leads  to  y  =  1.00  and  the  awful  result  x  =  0.00. 


equation  by  10  .  Then  the  system  would  be 


(14) 


"lO.O  x  +  100000  y  = 

_  1.00  x  +  1.00  y  = 
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We  see  that  poor  scal.ng  vitrc  a  wood  pivo’ nl  r.t i atcgy  forces  us  into  the 
same  enormous  rounding  error  thr  t  we  obtained  in  fee  8  from  the  original  set  of 
equations  and  a  bad  pivotal  strategy. 

The  conclusion  of  this  is  tha*  a  good  pivotal  strategy  is  only  good  when 
the  matrix  is  properly  scaled  in  ad v nee.  However,  it  must  be  admitted  that  so 
far  we  do  not  know  guaranteed  algorithms  for  scaling  matrices  well. 

It  is  normal  to  scale  matrices  by  simply  multiplying  rows  and  column  by 
factors.  In  effect,  one  chooses  nonsingular  diagonal  matrices  and 

nd  then  scales  A  by  the  t.ransforma’-ion 

A  -*  D’1  A 

Because  cond(A)  is  an  ingredient  of  all  our  error  bounds  and  convergence 
theorems,  it  is  natural  to  wish  to  select  and  so  as  to  reduce 

cond(D~1  A  D^)  to  as  low  a  value  as  is  reasonably  possible. 

One  usually  uses  powers  of  tne  floating-point  base  for  scnle  factors, 
to  avoid  the  introduction  of  rounding  errors  in  the  scaring.  Or,  alternatively, 
one  may  use  the  scaling  only  implicitly,  without  actually  altering  the  element. s 
of  A„ 

theorem  (F.  L,  Bauer)  If  the  ordered  set _ of  pivotal  elements  is  selected 

in  advance,  scaling  of  a  matrix  A  by  peters  of  the  floating-point  base  does  nc 
change  a  single  digit,  of  the  significand  of  any  intermediate  or  final  number 
in  the  solution  of  Ax  -  b  by  Gaussian  elimination. 

The  theorem  was  presented  in  Ba-er  [  1 .1  Thus  the  only  possible  effect 
of  the  scaling  of  A  on  the  rounding  errors  m^s’  occur  through  changing  the 
order  of  pivots.  Our  example  showed  that  the  change  in  pivots  can  indeed  make 
a  great  deal  of  difference. 

One  is  sometimes  advised  to  pick  and  so  that  the  resulting 
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matrix  A  has  its  maximum  element  m  each  row  and  column  (in  aoso'!  ..v- 

v^lue)  in  the  interval  [  ,1,  l),  in  whatever  number  bn  e  one  i".  using  .  However, 
Richard  Hamming  has  showed  (unpublished)  that  this  advice  does  not  always  lead 
to  good  scaling.  If 


A  = 


1 

2 


1 

-1 


2  x  10' 
109 


Then  both  of  the  following  matrices  are  decimally  scaled  equivalents  of  A: 


10-10 

o 

• — I 

1 

o 
* — 1 

2 

2  X  10' 10 

-io'10 

..1 

«  x 

^  2 

0 

However,  is  a  well-conditioned  matrix  that  offers  no  difficulties  :n  the 

solution  of  an  equation  sys-em,  whereas  A^  is  most  ill-conditioned  and 

provides  vast  troubles  for  elimination. 

Bauer  [2]  has  studied  the  problem  of  finding  2^  and  to  minimise 

cond(D*J'  A  2^).  It.  turns  out  that  the  solution  depends  on  cer'ain  properties 

of  the  nonnegative  matrices  1 A  |  •  !a  and  |  A  “*“ !  ■  |a|  .  (Here  j  B  | 

denotes  the  matrix  of  absolute  values  jb  .j.)  Clearly,  ve  can  nardly  hope  ‘o 

i  >  ,1 

compute  in  order  tc  find  a  reasonable  scaling,  so  that  we  can  compute  A 


oo,  it  is  an  open  question,  hov;  to  find  a  demons* r ably  good  and  convenient  scaling 
algorithm.  Existing  algorithms  are  either  very  superficial  or  potentially  very 
slow. 

The  only  cheerful  side  of  the  scaling  quest  ion  is  that  it  seems  to  be  a 
rare  matrix  which  good  scaling  changes  from  untractable  to  tractable  I 


11.  Analysis  of  rounding  errors.  We  pointed  out  in  Sec.  5  that  the 
direct  rounding  error  analysis  of  von  Neumann  and  Goldstine  was  extremely 
T  edious  to  apply.  Givens  [17J  introduced  the  idea  of  inverse  rounding  errors., 
Wilkinson  has  developed  this  into  a  very  powerful  tool  for  bounding  the  rounding 
errors  in  matrix  computations.  The  error  bounds  of  Secs.  8  and  9  were  obtained 
from  inverse  analysis.  Th^  basic  idea  is  to  change  the  nonassociative,  non- 
dist.ributive  floating-point  arithmetic  system  into  an  associative,  distributive 
number  system,  by  throwing  the  errors  back  onto  the  data  of  the  computation. 

For  example,  let  fl(u  X  v)  stand  for  the  floating-point  product 
^number  base  P  )  of  the  floating-point  numbers  u  and  v.  The  direct  error 
analysis  uses  statements  of  the  form 

w  =  fl(u  X  v)  --  uv  +  where  |:|  <  -  |uv|£I_t 

Farther  operations  on  w  introduce  new  errors,  and  one  has  to  keep  account  of 
the  cumulation  of  all  *he  old  and  new  rounding  errors.  Eventually.,  one  bounds  T  he 
difference  between  the  computed  final  answer  and  the  mathematically  correct 
i  nswer  corresponding  to  the  given  data. 

In  inverse  analysis,  one  makes  statements  of  the  form 

w  =  fl(u  X  v)  -  uv(l  +  6),  where  |&  |  <  |  p1”t  . 
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Thus  the  computed  product  is  considered  the  true  mathematical  product  of  (for 
example)  the  real  numbers  u  and  v(l  +  5),  which  differ  slightly  from  u 
and  v.  Further  floating-point  operations  on  w  produce  numbers  which  are 
always  treated  as  the  results  of  exact  operations  on  other  slightly  more 
perturbed  approximations  to  the  original  data.  The  final  answer  is  considered 
as  the  exact  solution  of  an  original  problem  with  data  which  are  perturbed  by 
amounts  for  which  bounds  are  given. 

If  desired,  these  inverse  error  bounds  can  be  converted  to  ordinary  erroi 
bounds,  by  normal  mathematical  methods. 

Inverse  error  analysis  turns  out  to  be  extremely  well  adapted  to  the 
analysis  of  algorithms  of  a  marching  type  which  continally  introduce  new  data. 

Both  the  solution  of  linear  equations  and  the  evaluation  of  polynomials  are  of 
this  type.  Inverse  error  analysis  is  not  at  all  well  suited  for  problems  of  an 
iterative  nature--for  example,  the  Newton  process  for  evaluating  the  snuare  root  of 
a  number. 

The  reader  is  referred  to  Wilkinson  [60,  62]  for  further  study  of  inverse 
round-off  analysis. 

A  second  .approach  to  round-off  analysis  is  the  interval  analysis,  ex¬ 
tensively  developed  by  Moore  [ h-o ] ,  but  based  on  the  idea  of  "range  numbers" 
presented  earlier  by  Dwyer  [ 6 J .  In  its  original  form,  interval  analysis  is 
poorly  adapted  to  matrix  computations,  but  Hansen  [253  has  modified  it  ingeniously 


for  matrix  work. 


12  Eigenvalues  of  symmetric  matrices. 


Space  does  not  permit  as 


extensive  a  treatment  of  the  eigenvalue  problem  as  that  given  for  the  linear 
equations  problem.  We  can  only  mention  a  few  highlights  of  today's  methods 
The  reader  is  referred  to  Wilkinson's  treatise  [6l]  for  an  almost  complete 
presentation  of  the  state  of  the  art. 

As  with  the  linear  equations  problem,  the  computation  of  eigenvalues 
of  matrices  divides  into  two  methods,  according  to  the  nature  of  the  matrices. 
For  large,  sparse  matrices  the  methods  are  mostly  infinite  iterations,  and  will 
not  be  considered  here.  For  dense,  stored  matrices,  most  methods  are  finite 
algorithms . 

If  a  matrix  A  is  symmetric,  its  eigenvalues  are  very  well  determined 
by  the  data.  In  fact,  let  the  symmetric  matrix  B  -  A  +  E  have  eigenvalues 

and  let  A  (also  symmetric)  have  eigenvalues  0^.  Then  the  eigenvalues  can 
be  so  numbered  that 

(15)  I  Oh  -  Pj  I  <  |lE||  ,  for  ail  1. 


Now  inverse  error  analysis  refers  the  computed  eigenvalues  of  a  matrix  A  back 
to  a  matrix  B  -  A  +  E.  If  E  can  be  proved  to  be  small  (as  it  can),  then  (l5..' 
shows  how  small  the  eigenvalue  errors  are.  In  fact,  today's  methods  can  yield 
eigenvalues  that  are  in  error  by  only  a  few  digits  in  the  least  significant  digits 
of  the  large  eigenvalues. 

The  method  of  Jacobi  [5l]  is  an  infinite  iteration  for  dense,  stored 
matrices.  It  produces  a  sequence  of  matrices  orthogonally  congruent  to  At 


A 


k 


Moreover, 


\ 


converges  to  a  diagonal  matrix 


D  whose  diagonal  entries  are,  of 
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course,  the  eigenvalues  of  A.  In  fact,  each  A. 


k+1 


is  computed  from  tne 


previous  by  a  rotation  in  the  coordinate  2-space  of  some  two  indices 

i  and  j,  a  rotation  chosen  so  that  =  0. 

J 

For  any  k  such  that  A^  is  almost  diagonal,  the  columns  of  the 
corresponding  orthogonal  matrix  U  are  approximately  column  eigenvectors  of  A 

K 

Moreover,  the  columns  are  themselves  orthogonal.  Thus  the  Jacobi  method  yields 
approximate  eigenvectors  of  fine  quality  as  a  by-product  of  the  basic  iter'Mon. 
The  whole  program  is  easy  to  write,  and  it  is  difficult  for  it  to  be  done  badly. 
There  are  some  theoretical  problems  about  how  good  the  eigenvectors  are,  and 
whether  the  actually  converge., 

Goldstine,  Murray,  and  von  Neumann  [19]  analyzed  the  rounding  errors  in 

a  fixed-point  version  of  the  Jacobi  method. 

The  original  Jacobi  algorithm  chose  i  and  j  to  maximize  the  absolute 
(k) 

value  of  the  element  a.  '  of  A,  .  Modern  algorithms  modify  this  criterion 

i>  J  k 

in  one  of  two  ways: 

(i)  In  the  cyclic  Jacobi  methods,  the  off-diagonal  elements  a,  are 

1>  J 

zeroed  in  some  cyclic  order.  Forsythe  and  Henrici  [ill  proved  the  convergence  o: 

a  common  cyclic  method.  See  also  Hansen  [22], 

(li)  In  threshold  Jacobi  methods,  an  element  a.  .  is  selected  for 

: 

annihilation  only  when  its  absolute  value  is  above  a.  certain  threshold  size, 
which  gets  smaller  as  the  iteration  progresses,  See  Fope  "nd  Tompkins  [k9J  1  r.d 
Corneil  l^]. 

It  has  been  proved  only  recently  that  the  cyclic  Jacobi  method  converges 


:  isdratically  for  any  matrix  A.  See  Schonhage  [52]  and  Wilkinson  [59].  "he 
work  was  based  on  that  of  Henrici  [2^]. 

Givens  [ 17 1  observed  that,  although  it  takes  an  infinite  sequence  of 


rotations  to  bring  to  diagonal  form,  a  mere  ~(n-l)(n-2)  rotations  can 

bring  to  tridiagonal  form.  This  reduced  the  problem  to  that  of  finding 
eigenvalues  of  tridiagonal  matrices,  and  the  latter  problem  has  been  a  subject 
of  research  ever  since.  See  Ortega  [43 J,  Ortega  -nd  Kaiser  [4-3],  and  recent  work 
of  Kahan  and  Varah  [33].  In  any  case,  the  Givens  idea  cut  the  practical  time 
of  finding  eigenvalues  by  a  factor  of  about  9  in  practice  (Wilkinson  [6l,  p.  335' 
A  few  years  later  Householder  (see  Householder  and  Bauer  [29])  introduced  a 
new  method  of  tridiagonalizing  a  symmetric  matrix,  using  n-2  reflections  instead 
of  ^(n-l)(n-2)  rotations.  This  cut  the  time  down  by  another  factor  of  two, 
and  effectively  put  the  Givens  method  out  of  business.  An  error  analysis  is 
given  by  Ortega  [44].  Most  contemporary  programs  use  the  Householder  method, 
but  differ  widely  in  how  eigenvalues  of  tridiagonal  matrices  are  found.  Getting 
the  eigenvectors  is  surprisingly  tricky,  and  lack  of  knowledge  of  how  to  do  it 
is  one  reason  for  the  occasional  continued  use  of  the  Jacobi  methods. 

13.  Eigenvalues  of  unsymmetric  matrices.  The  area  of  greatest  activity 
in  the  past  decade  of  research  on  computational  linear  algebra  hfs  been  the 
eigenvalue  problem  for  unsymmetric  matrices.  Only  one  method  from  before  the 
computer  era  is  still  in  use--the  power  method--and  it  has  only  limited  appli¬ 
cations  today.  Most  methods  in  use  today  were  unheard  of  15  years  ago. 

It  is  essential  to  realize  tne  instability  inherent  in  the  eigenvalue 
problem  for  unsymmetric  matrices.  In  contrast  to  the  close  bound  (l5),  for 
unsymmetric  matrices  the  corresponding  result,  due  to  Ostrowski  [46],  is 

06)  | on  -  Pi|  <  (polynomial  in  n)  X  ||e|| (all  i)  . 

The  above  result  is  very  weak,  and  yet  is  the  best  possible  general  result 
of  its  kind.  For  the  matrix 
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of  order  n  has  all  eigenvalues  0  for  e  =  0,  but  all  eigenvalues  •  d;st;r. 
h-ve  modulus  ld1/n  for  e  ^  0  Thus,  if  n  =  4o  and  e  =  10  r>j  .  eigenva 

have  modulus  0.1  (l). 

Fortunately,  eigenvalues  are  not  usually  so  sensitive  In  fact,  nifferen 
eigenvalues  of  a  matrix  A  c^n  differ  enormously  in  their  sensitivity  to 
perturbations  in  A.  Chapter  2  of  Wilkinson  [6l]  is  full  of  useful  results. 
They  are  generally  a  posteriori  results,  giving  bound  for  tne  changes  .r.  eigen¬ 
values  as  functions  of  perturbations  in  a  matrix  and  information  about  the 
other  eigenvalues  and  eigenvectors. 

The  great  power  of  the  Jacobi  method  for  symmetric  matrices,  '>nd  tne 
extremely  pleasant  rounding  characteristics  of  'unitary  matrices  led  to  n  d-.  sire 
to  use  them  for  the  unsymmetric  eigenvalue  problem.  Tne  b-  sic  tneore-r  is  d^<= 

• o  Schur : 

Theorem.  For  an  arbitrary  matrix  A,  r  here  exists  a.  unitary  nr.  '.ox  V 
such  that 

T  =  \l\  U 

is  triangular.  (Here  denotes  the  conjugate  transpose  of  U.) 

Since  the  eigenvalues  of  A  are  the  diagonal  elements  of  T,  the  hope 
has  been  to  find  unitary  matrices  which  bring  A  nearly  into  a  triangular  fcr^ , 
and  then  let  the  diagonal  elements  serve  as  approximate  eigenvalues  of  A 
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Investigations  by  Gree-nstadt  [2l]f  Lotkin  [ 37 J >  and  Eberlein  [ 7 J  offer 

some  hope,  but  no  real  promise  of  success 

For  most  methods  of  attacking  tne  eigenvalue  problem,  the  first  step 

is  to  condense  the  data,  to  save  time  and  storage  in  further  work  The  now 

universally  accepted  condensed  form  is  the  Hessenberg  matrix,  in  which  a,  .  =  0 

*  .  i ,  J 

for  i  -  j  >  1  (or  its  transpose)  It  is  possible  to  transform  A  by  orthogonal, 
congruences  of  the  Householder  type  into  -  Hessenberg  form  with  only  very 
small  rounding  errors.  Any  further  condensation  (say,  into  tridiagonal  form) 
is  subject  to  serious  losses  of  digits.  A  transformation  to  the  companion- 
matrix  form  is  particularly  disastrous  m  practice,  and  it  noriiiaily  requires  very 
substantial  increases  in  precision  to  successfully  yield  the  eigenvalues  of  A 
As  an  alternative,  one  can  transform  A  to  Hessenberg  form  by  Gaussian 
elimination  with  partial  pivoting,  a  similarity  transformation 

The  next  stage  in  the  unsymmet ric  eigenvalue  problem  is  to  get  the 
eigenvalues  of  a  Hessenberg  matrix  H.  A  v°nety  of  methods  h' ve  been  used. 

(i)  One  can  search  for  zeros  of  det (H  -  z! )  by  root-finding  methods, 
for  complex  z~  The  most  satisfactory  method  appears  to  be  that,  proposed  by 
Hyman  [30,  and  developed  by  Par  let  t  [1+7.  in  a  nwroer  of  programs  He  makes 
use  of  the  method  of  Laguerre  to  find  the  zeros  of  f(z),  following  by  a  form  of 
zero  suppression  Very  satisfactory  recurrences  are  used  to  evaluate  f(z), 
f'(z),  and  f " ( z ) ,  as  needed  by  the  Lag. err e  process.  After  the  eigenvalues 
have  been  fcund,  they  are  suppressed  by  applying  the  Laguerre 

process  to 

f ( z )/  n  -  v  ■ 

i  =  1 

(ii)  The  LR-algorithm  of  Rutishauser  [58'  was  an  important  development, 

Since  it  has  now  been  pretty  weil  supplanted  by  the  QR- algorithm ,  we  emit 


mention  of  it . 
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(xii)  Francis  [l5,  161  in  England,  and  Kublanovskaja.  [35]  in  the  Soviet 
Union  devised  the  very  interesting  QR-algorithm.  This  is  now  widely  considered 
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the  most  satisfactory  eigenvalue  algorithm  for  dense,  stored  unsymmetrie  matrices 

The  basic  theorem  is  that  an  arbitrary  real  square  matrix  A  can  be 

factored  in  the  form  A  =  QR,  where  Q  is  orthogonal,  and  where  K  is  an 

upper-triangular  matrix  with  all  diagonal  elements  r..  .  nonnegative.  If  A 

i> 1 

is  nonsingular,  then  both  Q  and  R  ^re  unique. 

In  f°ct,  the  computation  is  done  by  building  up  an  orthogonal  matrix 
T  T 

Q1  such  that  Q  A  =  R,  where  R  has  the  above  properties. 

As  an  aside,  for  nonsingular  A,  the  reader  will  be  more  familiar  with  the 

stepwise  determination  of  an  upper-triangular  matrix  R  with  positive  r^  ^  sech 

that  AR is  an  orthogonal  matrix  Q.  This  is  the  matrix  expression  of  the 

familiar  Gram-Schmidt,  process  of  analysis.  It  will  perhaps  surprise  the  reader 

that  the  matrix  Q  resulting  from  the  Gram-Schmidt  algorithm  is  normally  far 

from  orthogonal,  because  of  rounding  errors.  On  the  other  hand,  if  the  same  0 

T 

is  determined  so  that  Q  A  =  R,  the  rounding  errors  are  very  small. 

The  basic  QR-algorithm  proceeds  as  follows.  Let  H  =>  be  a  Hesser.be rg 

matrix-  For  k  =  0,  1,  2,  ...,  factor  in  the  form 


and  then  form 


■  RA 


It  is  easily  shown  that  H  is  also  a  Hessenberg  rr.aorix. 
The  basic  theorem  is  the  following. 

Tneorem,  Let  H  have  eigenvalues  A.^,  A.^,  .... 


with 

(17) 


\\  |  <  |\  |  <  ...  <  |\  | 

1  d  n 
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Then  the  matrices  converge  to  an  upper-trr  ngular  matrix  whose  diagonal 

—  —  K  ■  ■  -  — i  i  r|1  -  -  —  . r 

elements  are  the  eigenvalues  of  H 

In  the  more  usual  case  where  (17)  is  not  satisfied,  we  find  that 

converges  in  shape  to  a  blockwise  triangular  matrix,  (This  means  that  outside 

a  blockwise  triangular  matrix  all  elements  of  H  tend  to  zero,  as  k  -*  00  , 

but  ttr t  some  elements  of  the  blockwise  triangular  form  may  not  converge.) 

Moreover  the  2-by-2  and  1-by-l  diagonal  blocks  of  the  matrix  have 

eigenvalues  which  in  their  totality  converge  to  the  eigenvalues  of  H„ 

For  simplicity,  consider  a  matrix  H  with  eigenvalues 

0  <  V,  <  V,  <  •  ‘ 

12  n 

The  OR  method  for  such  a  matrix  converges  with  ran  error  which  is 


ot(VjA2)k] 

The  convergence  would  be  more  rapid  if,  instead  of  H,  we  dealt  with  the 
matrix  H  -  pi,  where  0  <  p  <  ^ .  If  p  were  practically  equal  to 
the  convergence  would  be  extremely  rapid.  Modifications  of  the  Of -algorithm 
have  been  devised  that  simulate  this  so-called  origin  shift  which  introduces 
p  near  \  .  After  one  eigenvalue  has  been  isolated,  the  QR  method  can  th 

be  applied  to  an  n-1  by  n-1  matrix  with  eigenvalues  . ..,  ,  New 

origin  shifts  are  then  introduced  to  bring  out  ^  as  rapidly  as  possible.  E‘ 
With  well  devised  origin  shifts,  the  whole  process  has  been  observed  to  convert 
with  an  average  of  less  than  two  iterative  steps  per  eigenvalue. 

Most  research  goes  into  the  invention  of  origin  shifts  when  some  of  the 
eigenvalues  are  complex  and  of  equal  modulus.  We  shall  not  attempt  to  give  i  .r 
ideas . 

A  more  recent  convergence  proof  has  been  given  by  Wilkinson  [62J,  but. 
like  the  Francis  proof,  is  given  for  an  arbitrary  matrix  A.  If  one  limits 
himself  to  matrices  of  Hessenberg  form,  easier  proofs  can  be  given;  see 
Kahan  (unpublished). 
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Normally,  the  eigenvalues  are  obtained  in  order  of  increasing  modulus., 

^arlett  [^8]  has  given  theorems  stating  precisely  when  this  occurs. 

If  H  is  a  symmetric  band  matrix,  then  the  QR-algorithm  preserves 
the  band  width  during  the  iteration,  and  is  very  satisfactory.  In  particular, 

9R  is  a  possible  algorithm  for  computing  eigenvalues  of  a  symmetric  tridiagonnl 
matrix. 

If  H  is  an  unsymmetric  band  matrix,  the  QR- algorithm  loses  the  zero  bands, 
above  the  diagonal. 

So  far,  we  have  not  mentioned  getting  the  eigenvectors  of  a  Hessenoerg 

matrix  H.  This  is  the  most  difficult  problem  we  shall  mention.  The  prevailing 

method  is  that  of  inverse  iter at  ion  .  The  eigenvalues  are  assumed  already  known. 

For  any  fixed  eigenvalue  one  selects  a  vector  arbitrarily.,  Then  one 

carries  out  an  iteration  of  the  following  form: 

For  each  k  =  0,  1,  2,  ...,  find  x  by  solving  4  he  system 

•  K+_L 

(18)  (H  -  \  I)xk+1  =  ^  . 

One  continues  until  is  quite  large.  In  easy  cases,  x^  is  nearly  an 

eigenvector  belonging  to  Wilkinson  [6l,  Chap,  9 j  discusses  variants  of  tms 

process.  Varah  [961  has  written  several  algorithms 

If  H  is  a  real  Hessenberg  matrix,  but  ^  is  a  complex  e.g.er.  value,  one 
has  to  choose  between  doing  complex  arithmetic,  or  ssrr.e  judicious _y  selected 
process  with  real  arithmetic. 

Finally,  one  transforms  x  back  to  the  original  coordinate  system  c.f  A 

K 

by  undoing  the  orthogonal  transformat  ions  from  A  to  H. 

If  some  of  the  eigenvalues  of  H  are  very  close,  the  real  problems  begin. 

A  pair  of  close  eigenvalues  may  in  fortunate  cases  have  distinct  column  eigenvecto- 
that  are  far  from  parallel;  this  represents  an  approximation  tea  double  eigenvalue 
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with  a  linear  elementary  divisor.  It  is  far  more  likely  that  a  pair  of  close 
eigenvalues  will  have  column  eigenvectors  that  are  almost  parallel  This  re¬ 
presents  an  anproximation  to  the  infinitely  more  probable  case  of  a  double 
eigenvalue  witn  nonlinear  divisors. 

In  the  former  case,  it  is  not  difficult  to  compute  two  eigenvectors  that 
are  far  from  parallel.  It  is  only  necessary  to  carry  out  the  iteration  (l8) 
with  different  x^,  or  with  two  slightly  different  values  of 

In  the  latter  case,  it  appears  difficult  to  obtain  much  from  the  iteration 
but  a  single  eigenvector  belonging  to  What  should  be  done  next?  In  parr 

one  doesn't  know  what  the  problem  proposer  would  like.  In  part,  one  doesn't 
know  what  is  possible.  Varah  is  carrying  out  research  on  the  problem.  He  is 
attempting  to  find  an  orthogonal  basis  for  the  invariant  subspace  of  dimension  2 
(in  this  case)  belonging  to 

For  a  "nice"  matrix,  Varah  is  also  getting  guaranteed  error  bounds  for 
all  eigenvalues  and  all  eigenvectors,  using  Gerschgorin  theorems,  as  Wilkinson 
recommends . 

l4.  Conclusion  and  moral..  The  computational  methods  of  linear  algebra 
are  moving  into  a  stage  where  we  have  reasonably  satisfactory  methods  for  dense, 
stored  matrices  A„  The  main  exception  is  the  problem  of  getting  eigenvectors 
with  error  bounds,  for  unsymmetnc  matrices.  The  algorithms  have  been  refined 
several  times,  and  are  being  published,  particularly  in  Numerisc he  Mat hematiK . 
Casual  users  of  matrix  algebra  will  do  no  better  thr,n  to  take  such  algorithms 
"off  the  shelf"  for  their  problems.  The  best  algorithms  are  mainly  written  in 
Algol  60.  Even  though  the  reader  may  use  another  language,  it  is  unquestionably 
worthwhile  for  him  to  learn  to  read  Algol  60,  just  in  order  to  be  able  to  read 
these  algorithms  and  adapt  them  to  his  own  problems. 
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No  method  of  solving  n  computet  1  on?  1  problem  is  re'  lly  rri'sb  -  a 
user  until  it  is  completely  described  in  an  algebraic  comput  ing  language  and 
made  completely  reliable.  Before  that,  there  are  indeterminate  aspects  :n  every 
algorithm,  Freauently  the  entire  advantage  of  a  certain  computing  process  .  les 
:n  t.he  treatment  of  certain  fine  points  which  can  hardly  be  suspec’ed  an4  :  i  :  ney 
are  completely  programmed.  This  is  the  reason  why  the  amateur  should  either 
consult  an  expert,  or  take  great  pains  to  pick  up  a  foolproof  'lgoritnm  'o.s 
is  the  reason  why  professionals  should  concentrate  ver y  hard  on  comp  let  e  .y 
foolproofing  the  algorithms  they  devise,  before  pitting  them  on  the  c.-.elf  : u • 
widespread  use. 
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